#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)

source("cueing_functions.R")

## This setup allows the function to be called from the
#  command line, making it easy to experiment with
## other thresholds for cocs and oppcs.
cocs <- as.numeric(args[1])
oppcs <- as.numeric(args[2])

# Set up data for use with ASA standard errors
pt_x <- read.csv(paste0("pt_cueing_data1616splitweighted.csv"))
pt_x$deltay <- (pt_x$agree2 - pt_x$agree1)
pt_x$friend <- apply(cbind(pt_x$copartisans, pt_x$cs), 1, 
                  function(z) ifelse(z[1] == 1, ifelse(z[2] >= cocs, 1, 0), 
                                     ifelse(z[2] >= oppcs, 1, 0)))

pt_x$legA2 <- apply(pt_x, 1, function(z) paste0(z[16], z[7]))
pt_x$legB2 <- apply(pt_x, 1, function(z) paste0(z[17], z[7]))

Aind <- which(colnames(pt_x) == "legA2")
Bind <- which(colnames(pt_x) == "legB2")

pt_x$dyads <- apply(pt_x, 1, function(r) paste0(r[Aind], "-", r[Bind]))

# Model is here 
fmla <- formula(deltay ~ treat + friend +
                  I(treat * copartisans) + 
                  I(treat * friend) + 
                  I(friend * copartisans) + 
                  #I(sgcs + ogcs) + I(sfcs + ofcs) + 
                  #I(treat * (sgcs + ogcs)) +
                  #I(treat * (sfcs + ofcs)) + 
                  sgcs + I(treat * sgcs) +
                  ogcs + I(treat * ogcs) + 
                  sfcs + I(treat * sfcs) + 
                  ofcs + I(treat * ofcs) +
                  interaction(factor(congress), committee,
                              leg_party, copartisans) + 
                  n1 + n2)

pt_fit <- lm(fmla, data = pt_x, weights = pt_x$weights)
coef(pt_fit)[1:14]
coef(pt_fit)[1:14]/sqrt(diag(vcov(pt_fit))[1:14])

pt_index <- unique(c(as.character(pt_x$legA2), as.character(pt_x$legB2)))
pt_dyad.mat <- apply(pt_x[,c(Aind, Bind)], 2, as.character)

## You will need access to a cluster computer to calculate the standard errors
cl <- makeCluster(32)
registerDoSNOW(cl)
pt_dyad.vcov <- dyad.robust.se(pt_x, pt_fit, pt_index, pt_dyad.mat)
stopCluster(cl)

# Store with qr = FALSE to save memory
pt_fit <- lm(fmla, data = pt_x, weights = pt_x$weights, qr = FALSE)
save.image(paste0("pt_cueingresults_", cocs, oppcs, ".RData"))
